Libraries:
# libraries -------------------
library(dplyr)
library(tidyverse)
library(ggplot2)
library(GGally)
library(glmmTMB)
library(TMB)
library(emmeans)
library(DHARMa)
library(lmtest)
library(ggeffects)
library(performance)
Set working directory
knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')
Source sapling data:
source("Scripts/SA_Import.R")
Pivot wider to create dataframe where each row is for one plot and has total details for each species group
sa_merge2 <- sapling_merge %>%
pivot_wider(names_from = Species_Groups, values_from = Total_Tally)
Import time since data and add it to the PIRI dataset
time_since <- read_csv("CleanData/Treat_Year_Data.csv")
sa_merge3 <- merge(sa_merge2, time_since, by = 'Site')
#log transform time from treatment data
sa_merge3$l.TFT <- log(sa_merge3$Time_from_Treat)
Run the ‘Add_BA’ script and merge with dataset:
source("Scripts/Add_BA.R")
# merge with sapling dataset -------------------
sa_merge4 <- merge(sa_merge3, prism_BA, by = 'Plot_No')
Run ‘Ground_Data.R’ script and add it to sapling dataset:
source("Scripts/Ground_Data.R")
# merge with sapling dataset -------------------
sa_merge5 <- merge(sa_merge4, ground3, by = 'Plot_No')
Source and add veg data:
source("Scripts/Veg_Data.R")
# merge with sapling dataset -------------------
sa.all <- merge(sa_merge5, veg3, by = "Plot_No")
rm(sa_merge3,
sa_merge4,
sa_merge5)
Sapling data will be convered into 1m2 plots in order to compare across and reduce the amount of scales of data collection, reducing it to two: 1m2 plots and per hectare observations
Create log transformed categories for newly added variables, then select for just the desired variables:
sa.all$PIRI.1m <- sa.all$PIRI/25
sa.all$SO.1m <- sa.all$Shrub_Oak/25
sa.all$Other.1m <- sa.all$Other/25
sa.all$l.PIRI1 <- log(sa.all$PIRI.1m + 1)
sa.all$l.SO1 <- log(sa.all$SO.1m + 1)
sa.all$l.other1 <- log(sa.all$Other.1m + 1)
sa.all2 <- sa.all %>%
select(Treat_Type, Region, Site, Plot_No, PIRI, l.PIRI1, Shrub_Oak, l.SO1, Other, l.other1, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>%
arrange(Treat_Type)
Create a dataframe keeping sapling at 25m2 observation levels & log transforming them
sa.all3 <- sa.all
sa.all3$l.PIRI <- log(sa.all3$PIRI + 1)
sa.all3$l.SO <- log(sa.all3$Shrub_Oak + 1)
sa.all3$l.other <- log(sa.all3$Other + 1)
sa.all3 <- sa.all3 %>%
select(Treat_Type, Region, Site, Plot_No, PIRI, l.PIRI, Shrub_Oak, l.SO, Other, l.other, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>%
arrange(Treat_Type)
Select just for numerical vs log and then look at paired plots for sapling transformed to 1m2 observation level:
#not log transformed (1m2)
sa.num <- sa.all2 %>%
select(PIRI.1m, S0.1m, Other.1m, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)
ggpairs(sa.num)
ggpairs(sa.num, aes(color = Treat_Type))
#log transformed (1m2)
sa.numl <- sa.all2 %>%
select(l.PIRI1, l.SO1, l.other1, l.TFT, l.BA_HA, l.BA_piri, l.Mineral, avgLD_l, l.Veg_Total, Treat_Type)
ggpairs(sa.numl)
ggpairs(sa.numl, aes(color = Treat_Type))
rm(sa.num,
sa.numl)
Can see the correlation coefficients for linear (Pearsons) relationships. None of them appear very strong, except for ones that are analogs (avg LD vs mineral soil exposure; ba/ha vs piri ba/ha)
No significant correlations obvious from plots
Looking at pairwise again but with sa.all3 dataset (sapling counts at 25m2 observation level)
#not log transformed (25m2)
sa3.num <- sa.all3 %>%
select(PIRI, Shrub_Oak, Other, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)
ggpairs(sa3.num)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggpairs(sa3.num, aes(color = Treat_Type))
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#log transformed (25m2)
sa3.numl <- sa.all3 %>%
select(l.PIRI, l.SO, l.other, l.TFT, l.BA_HA, l.BA_piri, l.Mineral, avgLD_l, l.Veg_Total, Treat_Type)
ggpairs(sa3.numl)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggpairs(sa3.numl, aes(color = Treat_Type))
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
rm(sa3.num,
sa3.numl)
Again no strong correlations
Going to try some scatterplots:
#PIRI & SO
ggplot(sa.all3, aes(x=l.PIRI, y = l.SO))+
geom_point()+
geom_smooth(method = lm)+
facet_grid(~Treat_Type)
## `geom_smooth()` using formula = 'y ~ x'
# negative for control, positive for fallrx, flat for other
#PIRI and other
ggplot(sa.all3, aes(x=l.PIRI, y = l.other, fill = Treat_Type))+
geom_point()+
geom_smooth(method = lm) +
facet_grid(~Treat_Type)
## `geom_smooth()` using formula = 'y ~ x'
# negative for control & fall rx, weak negative for springrx and mow, positive for harvest
#PIRI and veg
ggplot(sa.all3, aes(x=l.PIRI, y = l.Veg_Total, fill = Treat_Type))+
geom_point()+
geom_smooth(method = lm) +
facet_grid(~Treat_Type)
## `geom_smooth()` using formula = 'y ~ x'
# positive for harvest, weak negative for springrx, negative for mowrx, fallrx, and control
#relationships varying by treatment type may be a reason why the model are having a hard time in DHARMa
major personal breakthrough for dummies. For zero inflation, lets see if any TT or region is lacking alltogether in PIRI or SO saplings
sa.all4 <- sa.all3 %>%
group_by(Region)
tapply(sa.all4$PIRI, sa.all4$Region, summary)
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.02459 0.00000 2.00000
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2013 0.0000 9.0000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.07273 0.00000 1.00000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.09574 0.00000 3.00000
# no PIRI in LI,
tapply(sa.all4$PIRI, sa.all4$Treat_Type, summary)
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2564 0.0000 9.0000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.04902 0.00000 2.00000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.05357 0.00000 2.00000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.05882 0.00000 2.00000
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01149 0.00000 1.00000
I have previously done model elimination with zero inflated poisson distro & with model w.o TT. Neither are better than zi nbinom2 models. zi ~1 has distribution issues, ~Region does not
Model has much better fit and zero inflation issues are resolved when treatment type variable is included. Time to do elimination again but with TT; all models fail with nbinom1 distro, but works with nbinom2 (these notes are from when I was using a zi poisson distro, using nbinom TT no longer seems to be important to model fit, but we will keep for the overarching story)
sa.m1b <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
#won't converge
# test piri ba
sa.m2 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m2) #248
check_collinearity(sa.m2) #not rank deficient
#test overall BA
sa.m2b <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m2b) #239.3 (AIC smaller when PIRI BA included instead of overall BA)
check_collinearity(sa.m2b) # this model is not rank deficient
# no ba to compare with ba/piri ba models
sa.m3 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m3) #252.2
check_collinearity(sa.m3) #not rank deficent
lrtest(sa.m2b, sa.m3) #p=0.0001, keep piri ba
lrtest(sa.m2, sa.m3) # p = 0.01
#I think these two variables will have issues of collinearity, but I would like to confirm that, but models with both won't converge or are rank deficient
#test SO
sa.m4 <- glmmTMB(PIRI ~ Treat_Type + l.other + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m4) #240
check_collinearity(sa.m4) #rank deficient
emmeans(sa.m4, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
# ok, so rank deficiency shows up in pairwise as NaN
lrtest(sa.m2b, sa.m4) # p = 0.08 # drop
#test other
sa.m5 <- glmmTMB(PIRI ~ Treat_Type + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m5) #240.6
check_collinearity(sa.m5) #not rank deficient
#test mineral
sa.m6 <- glmmTMB(PIRI ~ Treat_Type + l.BA_piri + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m6) #238.8
check_collinearity(sa.m6) #not rank deficient
lrtest(sa.m5, sa.m6) # both other and mineral non-significant
#test veg
sa.m7 <- glmmTMB(PIRI ~ Treat_Type + l.BA_piri + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m7) #237.8
check_collinearity(sa.m7) #RANK DEFICIENT
lrtest(sa.m6, sa.m7) #drop
sa.m7b <- glmmTMB(PIRI ~ Treat_Type + l.BA_HA + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m7b) #247.7
check_collinearity(sa.m7b) # not rank deficient
#why is this happening
#if both BA/HA and piri BA/HA are important, should I make a proportion category? no - model using proportion won't converge
# sa.all3$BA_prop <- sa.all3$PIRI.BA_HA/sa.all3$BA_HA
# sa.all3$l.BA_prop <- log(sa.all3$BA_prop+1)
#
# hist(sa.all3$BA_prop)
#
# sa.m7p <- glmmTMB(PIRI ~ Treat_Type + l.BA_prop + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
# data = sa.all3,
# ziformula = ~Region,
# family = nbinom2)
#ok this idea failed. can't get models to work with piri ba proportion or log transformed piri ba prop
# i don't know how to resolve this issue or what the right steps to take
sa.m7c <- glmmTMB(PIRI ~ Treat_Type + l.BA_piri + avgLD_l + l.SO + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m7c) #237.3
# won't converge with l.other; above model is rank deficient
rm(sa.m7b, sa.m7c)
#test avgld
sa.m8 <- glmmTMB(PIRI ~ Treat_Type + l.BA_piri + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m8) #won't converge
check_collinearity(sa.m8) #RANK DEFICIENT
#test TT
sa.m9 <- glmmTMB(PIRI ~ l.BA_piri + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m9) #238.2
lrtest(sa.m9, sa.m2b) # p = 0.06
#test piri ba
sa.m10 <- glmmTMB(PIRI ~ l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m10) #245.2
#will avgld and piri ba converge without TT? NO
# tried with zi ~1, major dispersion issues when checking model fit
rm(sa.m1b, sa.m2, sa.m2b, sa.m3, sa.m4, sa.m5, sa.m10, sa.m11)
summary(sa.all3)
## Treat_Type Region Site Plot_No
## Length:498 Length:498 Length:498 Min. : 70.0
## Class :character Class :character Class :character 1st Qu.: 359.5
## Mode :character Mode :character Mode :character Median : 614.0
## Mean : 613.5
## 3rd Qu.: 887.5
## Max. :1143.0
## PIRI l.PIRI Shrub_Oak l.SO
## Min. :0.00000 Min. :0.00000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.0000
## Median :0.00000 Median :0.00000 Median :0.0000 Median :0.0000
## Mean :0.09438 Mean :0.04627 Mean :0.2691 Mean :0.1203
## 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:0.0000 3rd Qu.:0.0000
## Max. :9.00000 Max. :2.30259 Max. :7.0000 Max. :2.0794
## Other l.other Time_from_Treat l.TFT
## Min. :0.0000 Min. :0.000 Min. : 1.0 Min. :0.000
## 1st Qu.:0.0000 1st Qu.:0.000 1st Qu.: 3.0 1st Qu.:1.099
## Median :0.0000 Median :0.000 Median : 5.0 Median :1.609
## Mean :0.3514 Mean :0.173 Mean :10.9 Mean :1.804
## 3rd Qu.:0.0000 3rd Qu.:0.000 3rd Qu.: 8.0 3rd Qu.:2.079
## Max. :9.0000 Max. :2.303 Max. :33.0 Max. :3.497
## BA_HA l.BA_HA PIRI.BA_HA l.BA_piri
## Min. : 0.00 Min. :0.000 Min. : 0.0 Min. :0.000
## 1st Qu.: 7.00 1st Qu.:2.079 1st Qu.: 5.0 1st Qu.:1.792
## Median :16.00 Median :2.833 Median :11.0 Median :2.485
## Mean :16.47 Mean :2.521 Mean :12.6 Mean :2.218
## 3rd Qu.:23.00 3rd Qu.:3.178 3rd Qu.:18.0 3rd Qu.:2.944
## Max. :51.00 Max. :3.951 Max. :44.0 Max. :3.807
## Mineral_Soil l.Mineral Litter_Duff avgLD
## Min. : 0.0000 Min. :0.0000 Min. : 0.00 Min. : 0.000
## 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:90.50 1st Qu.: 2.250
## Median : 0.0000 Median :0.0000 Median :90.50 Median : 3.875
## Mean : 0.7199 Mean :0.2032 Mean :88.78 Mean : 4.179
## 3rd Qu.: 0.0000 3rd Qu.:0.0000 3rd Qu.:90.50 3rd Qu.: 5.500
## Max. :50.5000 Max. :3.9416 Max. :90.50 Max. :13.250
## avgLD_l Veg_Total l.Veg_Total
## Min. :0.000 Min. : 4.00 Min. :1.386
## 1st Qu.:1.179 1st Qu.: 37.50 1st Qu.:3.624
## Median :1.584 Median : 58.00 Median :4.060
## Mean :1.519 Mean : 62.64 Mean :3.970
## 3rd Qu.:1.872 3rd Qu.: 83.38 3rd Qu.:4.423
## Max. :2.657 Max. :202.00 Max. :5.308
library(caret)
## Loading required package: lattice
##
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
##
## lift
sa_check <- sa.all3 %>%
select(l.PIRI, l.other, l.TFT, l.BA_piri, l.Mineral, l.Veg_Total, avgLD_l)
findLinearCombos(sa.all3[,4:23]) #shows no linear combos. Don't get how the rank deficiency is happening
## $linearCombos
## list()
##
## $remove
## NULL
tapply(sa.all3$Other, sa.all3$Region, summary) #all present
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.3443 0.0000 6.0000
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 0.411 0.000 6.000
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2143 0.0000 6.0000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.6727 0.0000 9.0000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.3511 0.0000 4.0000
tapply(sa.all3$Shrub_Oak, sa.all3$Region, summary) #none in LI
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01639 0.00000 1.00000
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.000000 0.000000 0.006494 0.000000 1.000000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.000 1.964 3.000 7.000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2447 0.0000 6.0000
tapply(sa.all3$BA_HA, sa.all3$Region, summary) #all present
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 11.00 14.98 23.00 51.00
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.00 11.00 21.00 19.81 28.00 39.00
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 9.00 14.00 15.32 21.00 44.00
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 13.50 23.00 22.13 31.00 44.00
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 11.00 14.37 21.00 48.00
tapply(sa.all3$PIRI.BA_HA, sa.all3$Region, summary) #all present
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 2 9 10 16 41
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.00 11.00 16.00 17.41 23.00 37.00
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 9.00 11.44 16.00 34.00
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 8.00 18.00 19.18 30.00 44.00
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 9.00 10.32 16.00 32.00
tapply(sa.all3$Veg_Total, sa.all3$Region, summary) #all present
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.00 51.50 71.00 70.96 91.88 166.50
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.00 23.00 39.00 44.52 61.00 128.00
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.50 33.75 52.75 56.93 74.38 161.50
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 14.50 41.25 64.00 62.97 79.75 114.00
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.00 46.00 62.25 75.09 96.62 202.00
tapply(sa.all3$Mineral_Soil, sa.all3$Region, summary) #all present
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 1.586 0.500 30.500
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 0.589 0.500 15.500
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.6623 0.0000 50.5000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01818 0.00000 0.50000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2021 0.0000 8.0000
tapply(sa.all3$avgLD, sa.all3$Region, summary) #all present
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.000 1.750 1.977 2.750 7.250
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.75 4.00 5.25 5.61 6.75 12.50
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.250 3.562 5.000 5.347 6.750 13.250
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.750 2.500 3.500 3.686 4.500 9.750
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.750 2.812 4.000 4.298 5.438 12.250
tapply(sa.all3$Other, sa.all3$Treat_Type, summary) #all present
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.8889 1.0000 6.0000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.4118 0.0000 9.0000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.1429 0.0000 2.0000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.06618 0.00000 2.00000
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.1379 0.0000 2.0000
tapply(sa.all3$Shrub_Oak, sa.all3$Treat_Type, summary) #none in mowrx or springrx
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.7094 0.0000 7.0000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.4902 0.0000 6.0000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01786 0.00000 1.00000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
tapply(sa.all3$BA_HA, sa.all3$Treat_Type, summary) #all present
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 18.00 28.00 26.78 37.00 51.00
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 9.00 16.00 17.44 23.00 44.00
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 11.00 12.11 18.00 34.00
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2.000 7.000 9.059 14.000 37.000
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 10.00 14.00 15.85 21.00 41.00
tapply(sa.all3$PIRI.BA_HA, sa.all3$Treat_Type, summary) #all present
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 9.00 18.00 17.05 25.00 41.00
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 11.00 14.45 20.25 44.00
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 5.000 9.000 9.804 14.500 30.000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2.000 7.000 7.816 11.000 37.000
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 9.00 14.00 13.75 18.00 34.00
tapply(sa.all3$Veg_Total, sa.all3$Treat_Type, summary) #all present
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.00 26.50 46.00 46.85 61.00 115.50
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 14.5 46.0 68.0 72.7 91.0 202.0
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 15.00 36.12 55.00 55.21 68.62 132.00
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.00 37.00 66.75 70.47 94.38 196.50
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 16.00 49.50 61.50 64.65 81.00 131.00
tapply(sa.all3$Mineral_Soil, sa.all3$Treat_Type, summary) #all present
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.000000 0.000000 0.008547 0.000000 0.500000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 0.201 0.000 8.000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 1.652 0.000 50.500
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 1.199 0.500 30.500
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.9368 0.5000 15.5000
tapply(sa.all3$avgLD, sa.all3$Treat_Type, summary) #all present
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.500 3.250 4.750 5.024 6.750 12.500
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.750 3.250 4.375 4.529 5.750 9.750
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.250 3.688 4.775 5.291 6.750 11.000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.500 2.750 2.816 3.812 13.250
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.500 1.500 3.500 4.046 5.250 11.500
sa.m6 <- glmmTMB(PIRI ~ Treat_Type + l.BA_piri + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
summary(sa.m6) #238.8
## Family: nbinom2 ( log )
## Formula:
## PIRI ~ Treat_Type + l.BA_piri + avgLD_l + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Zero inflation: ~Region
## Data: sa.all3
##
## AIC BIC logLik deviance df.resid
## 238.8 306.2 -103.4 206.8 482
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 0.00000002537 0.0001593
## Site (Intercept) 0.58921621688 0.7676042
## Number of obs: 498, groups: Plot_No:Site, 498; Site, 47
##
## Dispersion parameter for nbinom2 family (): 0.551
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.4666 2.1272 -2.570 0.01017 *
## Treat_TypeFallRx -0.1423 0.8538 -0.167 0.86759
## Treat_TypeHarvest 1.6832 1.1041 1.524 0.12740
## Treat_TypeMowRx 1.6446 0.8920 1.844 0.06522 .
## Treat_TypeSpringRx -0.5944 1.3264 -0.448 0.65403
## l.BA_piri 1.3397 0.3983 3.363 0.00077 ***
## avgLD_l -0.7821 0.6039 -1.295 0.19533
## l.Veg_Total -0.3636 0.3773 -0.964 0.33516
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.9116 0.9406 2.033 0.0421 *
## RegionLI 21.4369 30910.2106 0.001 0.9994
## RegionMA -20.0122 13973.6041 -0.001 0.9989
## RegionME -1.7758 1.3951 -1.273 0.2031
## RegionNH -1.5524 1.1818 -1.314 0.1890
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sa.m6_sr <- simulateResiduals(sa.m6, n = 1000, plot = TRUE) #looks good
testResiduals(sa.m6_sr) #passes
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.046286, p-value = 0.2364
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.85775, p-value = 0.688
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00072289156626506 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.046286, p-value = 0.2364
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.85775, p-value = 0.688
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00072289156626506 )
## 0
testZeroInflation(sa.m6_sr) #passes
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99994, p-value = 1
## alternative hypothesis: two.sided
emmeans(sa.m6, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
## Treat_Type response SE df asymp.LCL asymp.UCL
## Control 0.0361 0.0286 Inf 0.00765 0.170
## FallRx 0.0313 0.0253 Inf 0.00641 0.153
## Harvest 0.1942 0.1760 Inf 0.03288 1.147
## MowRx 0.1869 0.1217 Inf 0.05215 0.670
## SpringRx 0.0199 0.0247 Inf 0.00175 0.227
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## Control / FallRx 1.153 0.984 Inf 1 0.167 0.9998
## Control / Harvest 0.186 0.205 Inf 1 -1.524 0.5465
## Control / MowRx 0.193 0.172 Inf 1 -1.844 0.3483
## Control / SpringRx 1.812 2.403 Inf 1 0.448 0.9917
## FallRx / Harvest 0.161 0.186 Inf 1 -1.579 0.5106
## FallRx / MowRx 0.167 0.163 Inf 1 -1.839 0.3510
## FallRx / SpringRx 1.572 2.132 Inf 1 0.333 0.9973
## Harvest / MowRx 1.039 1.138 Inf 1 0.035 1.0000
## Harvest / SpringRx 9.753 14.468 Inf 1 1.535 0.5393
## MowRx / SpringRx 9.384 12.813 Inf 1 1.640 0.4718
##
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log scale
None significantly different from eachother
sa.piri1 <- ggpredict(sa.m6, terms = c("avgLD_l", "Treat_Type"))
#control and fallrx are almost completely identical - if I wanted to nudge/jitter
# sa.piri2 <- sa.piri1 %>%
# filter(group == "FallRx")
#
# sa.piri2$predicted <- sa.piri2$predicted+0.01
#
# sa.piri3 <- sa.piri1 %>%
# filter(group != "FallRx")
#
# sa.piri3 <- rbind(sa.piri3, sa.piri2)
ggplot(sa.piri1, aes(x, predicted, color = group))+
geom_line(linewidth = 1, show.legend = FALSE) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type")+
theme(axis.text = element_text(size = 10), axis.title = element_text(size = 14))+
labs(x = "Average leaf litter depth (log transformed)",
y = NULL)
ggsave(filename = "SA_PIRI_avgld.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)
sa.piri4 <- ggpredict(sa.m6, terms = c("l.BA_piri", "Treat_Type"))
#again control and fall rx overlap - if I wanted to nudge/jitter
# sa.piri5 <- sa.piri4 %>%
# filter(group == "FallRx")
#
# sa.piri5$predicted <- sa.piri5$predicted+0.01
#
# sa.piri6 <- sa.piri4 %>%
# filter(group != "FallRx")
#
# sa.piri6 <- rbind(sa.piri6, sa.piri5)
ggplot(sa.piri4, aes(x, predicted, color = group))+
geom_line(linewidth = 1, show.legend = FALSE) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type")+
theme(axis.text = element_text(size = 10), axis.title = element_text(size = 14))+
labs(x = "Average pitch pine basal area per hectare (log transformed)",
y = "Pitch pine saplings \n (corrected for time)")
ggsave(filename = "SA_PIRI_BA.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)
sa.piri5 <- ggpredict(sa.m6, terms = c("l.Veg_Total", "Treat_Type"))
ggplot(sa.piri5, aes(x, predicted, color = group))+
geom_line(linewidth = 1, show.legend = FALSE) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type")+
theme(axis.text = element_text(size = 10), axis.title = element_text(size = 14))+
labs(x = "Average understory vegetation cover (log transformed)",
y = NULL)
ggsave(filename = "SA_PIRI_veg.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)
Trying to plot with jtools package
Start from the beginning knowing that both Region and Treat Type have 0 areas
tapply(sa.all4$Shrub_Oak, sa.all4$Region, summary)
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01639 0.00000 1.00000
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.000000 0.000000 0.006494 0.000000 1.000000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.000 1.964 3.000 7.000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2447 0.0000 6.0000
# no SO in LI,
tapply(sa.all4$Shrub_Oak, sa.all4$Treat_Type, summary)
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.7094 0.0000 7.0000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.4902 0.0000 6.0000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01786 0.00000 1.00000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
# no SO in MowRX or SpringRx
rm(sa.all4)
so.sa5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + avgLD_l + l.BA_piri + l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')
#won't converge
so.sa6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + avgLD_l + l.BA_piri + l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.sa6) #334.7 - this is the only way the model will converge, to remove Mineral
## [1] 334.7261
lrtest(so.sa5, so.sa6)
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + avgLD_l +
## l.BA_piri + l.BA_HA + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + avgLD_l + l.BA_piri +
## l.BA_HA + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 19
## 2 18 -149.36 -1
so.sa7 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + avgLD_l + l.BA_piri + l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.sa7) #333.3
## [1] 333.3224
so.sa8 <- glmmTMB(Shrub_Oak ~ Treat_Type + avgLD_l + l.BA_piri + l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.sa8) #331.3
## [1] 331.3264
so.sa9 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.Mineral + avgLD_l + l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.sa9) #326.3 (finally able to add l.Mineral back in)
## [1] 326.2949
lrtest(so.sa8, so.sa9) #mineral very significant
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + avgLD_l + l.BA_piri + l.BA_HA + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.Mineral + avgLD_l + l.BA_HA + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 16 -149.66
## 2 16 -147.15 0 5.0315 < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
so.sa10 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.Mineral + l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.sa10) #324.3
## [1] 324.3074
lrtest(so.sa9, so.sa10)
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.Mineral + avgLD_l + l.BA_HA + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.Mineral + l.BA_HA + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 16 -147.15
## 2 15 -147.15 -1 0.0125 0.9111
so.sa11 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')
AIC(so.sa11) #won't converge
## [1] 324.2177
so.sa10_sr <- simulateResiduals(so.sa10, n = 1000, plot = T) #passes
testResiduals(so.sa10_sr) #passes
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.034621, p-value = 0.5893
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.87777, p-value = 0.77
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.002008032
## sample estimates:
## outlier frequency (expected: 0.000542168674698795 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.034621, p-value = 0.5893
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.87777, p-value = 0.77
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.002008032
## sample estimates:
## outlier frequency (expected: 0.000542168674698795 )
## 0
testZeroInflation(so.sa10_sr) #passes
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99318, p-value = 0.654
## alternative hypothesis: two.sided
emmeans(so.sa10, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
## Treat_Type response SE df asymp.LCL asymp.UCL
## Control 0.0000155 0.104 Inf 0.000000000000000222 Inf
## FallRx 0.0000390 0.263 Inf 0.000000000000000222 Inf
## Harvest 0.0000818 0.552 Inf 0.000000000000000222 Inf
## MowRx 0.0000000 0.000 Inf 0.000000000000000222 Inf
## SpringRx 0.0000000 0.000 Inf 0.000000000000000222 Inf
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## Control / FallRx 0 0 Inf 1 -2.621 0.0665
## Control / Harvest 0 0 Inf 1 -1.415 0.6178
## Control / MowRx 2699054976 79659759310326 Inf 1 0.001 1.0000
## Control / SpringRx 37242947 314408982995 Inf 1 0.002 1.0000
## FallRx / Harvest 0 1 Inf 1 -0.624 0.9713
## FallRx / MowRx 6809895316 200986873784945 Inf 1 0.001 1.0000
## FallRx / SpringRx 93966433 793274786998 Inf 1 0.002 1.0000
## Harvest / MowRx 14295160591 421906579510505 Inf 1 0.001 1.0000
## Harvest / SpringRx 197251968 1665222450358 Inf 1 0.002 1.0000
## MowRx / SpringRx 0 424 Inf 1 0.000 1.0000
##
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log scale
so.p1 <- ggpredict(so.sa10, terms = c("l.BA_HA", "Treat_Type"))
ggplot(so.p1, aes(x, predicted, color = group))+
geom_line(linewidth = 1) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type")+
theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10), axis.title = element_text(size = 14), plot.title=element_text(hjust=0.5, size=20))+
labs(x = "Average basal area per hectare (log transformed)",
y = "Total Count of Shrub Oak",
title = "Predicted Counts of Shrub Oak Saplings \n >/= 2.5 cm & < 10 cm DBH")
#very small effect; springrx and mowrx both have 0 SO saplings
so.p2 <- ggpredict(so.sa10, terms = c("l.Mineral", "Treat_Type"))
ggplot(so.p2, aes(x, predicted, color = group))+
geom_line(linewidth = 1) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type")+
theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10), axis.title = element_text(size = 14), plot.title=element_text(hjust=0.5, size=20))+
labs(x ="Average exposed mineral soil (log transformed)",
y = "Total Count of Shrub Oak",
title = "Predicted Counts of Shrub Oak Saplings \n >/= 2.5 cm & < 10 cm DBH") +
xlim(0, 0.5)
## Warning: Removed 25 rows containing missing values or values outside the scale range
## (`geom_line()`).
this is a weird one, because exposed mineral soil clearly doesn’t go over a certain percent; the median is 0 and the mean 0.2, the IQR 0 - 0. Maybe it looks better at a smaller x scale.
summary(sa.all3$l.Mineral)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2032 0.0000 3.9416
summary(sa.all3$Mineral_Soil)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.7199 0.0000 50.5000
so.sa1 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + avgLD_l + l.BA_piri + l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(so.sa1) #333.8
## [1] 333.7745
# check_collinearity(so.sa1) #rank deficient
# test avgld
so.sa2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.BA_piri + l.BA_HA + l.PIRI + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(so.sa2) #331.9
## [1] 331.8841
check_collinearity(so.sa2) #not rank deficient
## # Check for Multicollinearity
##
## * conditional component:
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## Treat_Type 1.06 [1.01, 1.30] 1.03 0.95 [0.77, 0.99]
## l.BA_piri 2.38 [2.10, 2.74] 1.54 0.42 [0.37, 0.48]
## l.BA_HA 2.48 [2.18, 2.86] 1.58 0.40 [0.35, 0.46]
## l.PIRI 1.06 [1.01, 1.30] 1.03 0.95 [0.77, 0.99]
## l.other 1.23 [1.13, 1.39] 1.11 0.81 [0.72, 0.88]
## l.Mineral 1.00 [1.00, Inf] 1.00 1.00 [0.00, 1.00]
# test piri ba
so.sa3 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.BA_HA + l.PIRI + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(so.sa3) #329.9
## [1] 329.9233
lrtest(so.sa2, so.sa3) #drop
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.BA_piri + l.BA_HA + l.PIRI + l.other +
## l.Mineral + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.BA_HA + l.PIRI + l.other + l.Mineral +
## offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 18 -147.94
## 2 17 -147.96 -1 0.0393 0.8429
# test piri
so.sa4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.BA_HA + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(so.sa4) #328
## [1] 327.9643
lrtest(so.sa3, so.sa4)
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.BA_HA + l.PIRI + l.other + l.Mineral +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.BA_HA + l.other + l.Mineral + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 17 -147.96
## 2 16 -147.98 -1 0.041 0.8396
# anything dropped from model won't converge, higher aic with ~1, won't converge with ~TT
check_collinearity(so.sa4) #not rank deficient
## # Check for Multicollinearity
##
## * conditional component:
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## Treat_Type 1.02 [1.00, 4.31] 1.01 0.98 [0.23, 1.00]
## l.BA_HA 1.07 [1.02, 1.28] 1.04 0.93 [0.78, 0.98]
## l.other 1.07 [1.02, 1.28] 1.03 0.94 [0.78, 0.98]
## l.Mineral 1.00 [1.00, Inf] 1.00 1.00 [0.00, 1.00]
so.sa4_sr <- simulateResiduals(so.sa4, n = 1000, plot = TRUE) #looks good
testResiduals(so.sa4_sr) #passes
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.031283, p-value = 0.7144
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.8769, p-value = 0.93
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000542168674698795 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.031283, p-value = 0.7144
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.8769, p-value = 0.93
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000542168674698795 )
## 0
testZeroInflation(so.sa4_sr) #passes
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99681, p-value = 0.902
## alternative hypothesis: two.sided
emmeans(so.sa4, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
## Treat_Type response SE df asymp.LCL asymp.UCL
## Control 0.00000260 0.0935 Inf 0.000000000000000222 Inf
## FallRx 0.00000602 0.2162 Inf 0.000000000000000222 Inf
## Harvest 0.00000958 0.3439 Inf 0.000000000000000222 Inf
## MowRx 0.00000000 0.0000 Inf 0.000000000000000222 Inf
## SpringRx 0.00000000 0.0000 Inf 0.000000000000000222 Inf
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## Control / FallRx 0 0 Inf 1 -2.149 0.1996
## Control / Harvest 0 0 Inf 1 -1.004 0.8539
## Control / MowRx 7591553109 315497435464353 Inf 1 0.001 1.0000
## Control / SpringRx 42869891 323998497992 Inf 1 0.002 1.0000
## FallRx / Harvest 1 1 Inf 1 -0.353 0.9967
## FallRx / MowRx 17556161776 729616711771682 Inf 1 0.001 1.0000
## FallRx / SpringRx 99140550 749276197656 Inf 1 0.002 1.0000
## Harvest / MowRx 27922185717 1160418410274454 Inf 1 0.001 1.0000
## Harvest / SpringRx 157678021 1191685840742 Inf 1 0.002 1.0000
## MowRx / SpringRx 0 239 Inf 1 0.000 1.0000
##
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log scale
None significantly different. No SO saplings in SpringRx or MowRx - not sure if I should deal with that differently?
sa.so1 <- ggpredict(so.sa4, terms = c("l.BA_HA", "Treat_Type"))
# if I wanted to jitter/nudge lines
# sa.so2 <- sa.so1
# sa.so2 <- sa.so2 %>%
# filter(group == "SpringRx")
#
# sa.so2$predicted <- sa.so2$predicted+0.05
#
# sa.so3 <- sa.so1 %>%
# filter(group != "SpringRx")
#
# sa.so3 <- rbind(sa.so3, sa.so2)
ggplot(sa.so1, aes(predicted,x, color = group))+
geom_line(linewidth = 1) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
coord_flip()+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type")+
theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10), axis.title = element_text(size = 14), plot.title=element_text(hjust=0.5, size=20))+
labs(y = "Average basal area per hectare (log transformed)",
x = "Total Count of Shrub Oak",
title = "Predicted Counts of Shrub Oak Saplings \n >/= 2.5 cm & < 10 cm DBH")
# such a small effect
sa.so2 <- ggpredict(so.sa4, terms = c("l.Mineral", "Treat_Type"))
ggplot(sa.so2, aes(predicted,x, color = group))+
geom_line(linewidth = 1) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
coord_flip()+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type")+
theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10), axis.title = element_text(size = 14), plot.title=element_text(hjust=0.5, size=20))+
labs(y = "Average exposed mineral soil per 1m^2^ (log transformed)",
x = "Total Count of Shrub Oak",
title = "Predicted Counts of Shrub Oak Saplings \n >/= 2.5 cm & < 10 cm DBH")
#VERY WEIRD, likely because mineral soil is so limited
sa.so3 <- ggpredict(so.sa4, terms = c("l.other", "Treat_Type"))
ggplot(sa.so3, aes(predicted,x, color = group))+
geom_line(linewidth = 1) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
coord_flip()+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type")+
theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10), axis.title = element_text(size = 14), plot.title=element_text(hjust=0.5, size=20))+
labs(y = "Total counts of other sapling species (log transformed)",
x = "Total Count of Shrub Oak",
title = "Predicted Counts of Shrub Oak Saplings \n >/= 2.5 cm & < 10 cm DBH")
# so small